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ABSTRACT 

A  problem  of  gasification  and  heating  of  a  solid  due  to  the  action  of  an 
external  energy  source  is  discussed.  The  problem  involves  a  moving  boundary 
when  the  solid  gasifies.  At  parts  of  the  boundary  where  gasification  is 
taking  place,  a  model  problem  looks  very  much  like  the  one-phase  Stefan 
problem  with  an  energy  source  at  the  moving  boundary.  However,  any  gas 
produced  is  assumed  to  blow  away  immediately,  and  there  is  no  conduction  to 
the  solid  from  the  outside,  even  when  the  surface  temperature  of  the  solid  is 
below  tne  gasification  temperature.  Accordingly,  if  the  temperature  is 
extended  to  a  function  defined  over  all  space  by  setting  it  equal  to  the 
gasification  temperature  outside,  the  temperature  will  not  necessarily  be 
continuous  at  the  boundary,  and  instead  a  Neumann  condition  may  be  satisfied 
there.  Also,  no  resolidification  is  possible,  so  that  the  region  occupied  by 
the  solid  cannot  increase.  Thus,  one  has  the  possibility  of  a  situation  in 
which  the  boundary  may  alternately  move  and  be  stationary.  A  generalized 
formulation  of  the  problem  is  given,  a  numerical  algorithm  is  described,  and 
computational  results  are  presented. 
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SIGNIFICANCE  AND  EXPLANATION 


The  classical  Stefan  prohlem  models  the  evolution  of  the  ice-water 
interface  in  a  melting  process  as  it  moves  under  the  influence  of  heat 
conduction  in  the  two  phases.  The  gasification  problem  arises  when  one 
irradiates  a  solid  and  vaporizes  a  portion  of  it,  with  the  gas  blowing  away. 
The  gasification  problem  differs  from  the  Stefan  problem  in  that  the  front  can 
move  in  only  one  direction,  and  the  gas  cannot  resolidify. 

In  contrast  to  the  case  with  the  Stefan  problem,  in  the  gasification 
problem  the  solid  boundary  can  stop  moving  and  start  cooling  down.  In  that 
case  the  nature  of  the  boundary  conditions  to  be  satisfied  changes.  A  key 
goal  of  this  paper  is  to  devise  an  algorithm  which  will  automatically  solve 
one  or  the  other  type  of  boundary  value  problem  as  it  becomes  the  relevant 
one . 


The  paper  presents  an  algorithm  which  solves  a  variety  of  boundary  value 

problems  without  explicitly  locating  the  boundaries  in  question.  Also 

included  is  a  means  of  treating  the  deposition  of  energy  on  a  surface  without 

location  of  the  surface.  Numerical  results  are  obtained  and  compared  with  the 

qualitative  properties  of  solutions  previously  anticipated.  Accession  For 
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The  responsibility  for  the  wording  and  views  expressed  in  this  descriptive 
summary  lies  with  MRC,  and  not  with  the  author  of  this  report. 


NUMERICAL  SOLUTION  OF  A  GASIFICATION  PROBLEM 


Joel  C.  W.  Rogers 


1.  Introduction 

The  problem  we  consider  is  of  the  following  sort:  A  solid  material  is 
exposed  to  an  external  source  of  heat.  Heat  is  conducted  internally  through 
the  solid,  and  when  the  solid  reaches  a  critical  temperature,  it  gasifies. 

The  gas  thus  formed  blows  away  and  no  longer  interacts  with  the  solid.  In  the 
version  of  the  problem  given  here,  the  external  heat  source  is  radiant  energy, 
to  which  the  solid  is  opaque  and  to  which  the  medium  external  to  the  solid  is 
transparent.  However,  other  types  of  heat  sources  may  be  treated  as  well;  in 
particular,  heat  sources  of  a  frictional  nature  at  the  surface  of  the  solid 
may  be  considered,  in  which  case  the  problem  is  better  known  as  an  "ablation" 
problem. 

In  the  next  section  we  will  proceed  with  a  careful  description  of  the 
phenomena  we  expect  to  find  associated  with  the  gasification  process,  and  in 
particular  we  will  analyze  the  ways  in  which  this  problem  resembles  and 
differs  from  the  classical  Stefan  problem.  On  the  basis  of  this  analysis,  we 
will  begin  development  of  a  mathematical  model.  A  third  section  will  discuss 
more  specifically  a  time-discretized  version  of  the  model,  the  treatment  of 
boundary  conditions,  etc.  Following  this,  there  will  be  a  brief  description 
of  a  numerical  quadrature  of  the  model,  as  well  as  computational  results  for  a 
particular  problem  whose  initial  and  boundary  data  will  have  been  chosen  so  as 
to  bring  about  a  solution  exhibiting  the  phenomena  which  are  anticipated  in 
the  second  section.  A  final  section  addresses  the  question  of  the  nature  of 
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the  dependence  of  solutions  of  the  time-discretized  problem  on  the  initial  and 
boundary  data.  We  show  that  this  dependence  is  not  generally  monotone  in 
character.  Various  weakened  types  of  monotonicity  which  may  hold  are  only 
conjectured. 

We  should  make  it  clear  that  we  have  not  answered  the  critical  question 
» 

of  convergence  of  the  algorithm  presented  here,  nor  the  questions  regarding 
stability  and  regularity  for  the  approximate  solutions  generated  by  the 
algorithm.  However,  we  would  not  have  presented  this  paper  unless,  bolstered 
by  confidence  in  the  essential  correctness  of  our  analysis  of  the  salient 
features  of  the  gasification  problem  and  the  numerical  results  that  we  had 
obtained,  we  felt  that  the  convergence  of  the  algorithm  to  a  solution  of  the 
problem  could  be  proved. 

With  regard  to  the  mathematical  theory  of  solutions  of  nonlinear 
parabolic  equations,  we  place  the  importance  of  this  problem  in  the  following 
context.  The  theory  of  the  classical  Stefan  problem  we  consider  to  be  in 
fairly  good  shape,  in  respect  both  to  the  basic  questions  of  existence  and 
uniqueness  of  solutions  and  of  effective  computational  methods  for  solution. 
For  systems  of  degenerate  parabolic  equations  in  several  space  variables,  the 
situation  is  more  complicated:  The  proofs  of  monotonicity  and  stability  for 
solutions  of  the  Stefan  problem  do  not  go  over  to  the  general  case  of  systems, 
and  indeed  for  some  systems  modeling  phase  transitions  with  solute  diffusion, 
the  experimental  evidence  and  some  linear  perturbation  theories  point  to  a 
high  degree  of  instability  and  complexity  for  the  boundaries  between  different 
phases  [71.  The  difficulties  arising  in  the  treatment  of  parabolic  svstems  in 
several  space  variables  are  in  some  resnects  similar  to  those  encountered  in 
the  study  of  hyperbolic  systems  in  several  space  variables  [f].  The 
gasification  probLem  belongs  to  what  mav  e  the  simplest  type  of  nontrivial 
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system.  There  are  two  dependent  variables,  the  enthalpy  u  and  a  quantity 


X  representing  che  fraction  of  material  gasified  at  each  point,  but  in  fact 
in  many  cases  x  is  an  explicit  algebraic  function  of  u,  in  which  case  the 
problem  can  be  put  in  Stefan-like  form.  Departures  from  this  condition  are 
characterized  by  the  fact  that  x  is  constant  in  time.  Accordingly,  an 
analysis  of  the  mathematical  chanqes  wrought  by  enlargement  of  the  Stefan 
problem  to  the  gasification  problem  may  be  very  informative  with  respect  to 
the  phenomena  to  be  expected  when  one  deals  with  more  general  parabolic 
systems. 

From  the  point  of  view  of  numerical  analysis,  there  are  twc  aspects  of 
our  treatment  which  may  be  of  interest.  The  first  is  the  means  pe  use  to 
deposit  radiant  energy  on  a  moving,  generally  irregular,  surface  (the  gas- 
solid  interface)  without  explicitly  tracking  the  surface.  The  second  feature 
is  that,  in  certain  regimes  of  the  initial  and  boundary  data,  the  problem 
looks  like  a  parabolic  problem  with  fixed  boundary  and  Neumann  boundary  data. 
In  that  regime  our  algorithm  will  solve  the  problem  without  explicitly 
locating  the  boundary.  This  is  to  be  compared  with  limiting  regimes  in  which 
the  classical  Stefan  problem  reduces  to  a  parabolic  problem  with  fixed 
boundary  and  nirichlet  data.  Certain  algorithms  which  have  been  given  for  the 
Stefan  problem  will  solve  this  problem  without  explicitly  locatina  the 
boundary  [2,  1,  41. 
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2.  Governing  Assumptions 

Typically  in  the  gasification  problem,  we  will  assume  that  the  heat 
conduction  in  the  interior  of  the  solid  may  be  described  by  the  equation 

ufc  =  Af ( u )  ,  (2.1) 

where  we  may  think  of  u  as  the  "enthalpy"  per  unit  volume  and  f(u)  as  the 
"temperature".  f(u)  will  be  a  nondecreasing  Lipschitz-continuous  function 
of  u.  We  shall  assume  that  the  solid  material  will  only  gasify  upon  reaching 
a  critical  gasification  temperature,  which  we  may  take  equal  to  0.  Thus, 
(2.1)  will  hold  in  the  region  where  the  temperature  is  below  this  number,  that 
is,  is  negative. 

Suppose  the  energy  transferred  to  the  solid  by  the  external  energy  source 
is  F  times  a  Dirac  measure  on  the  surface.  When  the  temperature  of  the 
solid  material  at  the  surface  remains  below  the  gasification  temperature, 
energy  conservation  at  the  solid  boundary  is  expressed  as 

<f(u>)  =  F  ,  (2.2) 

dn 

where  n  is  the  unit  outward  normal  of  the  surface.  In  this  case  the 
boundary  is  not  changinq,  since  no  gasification  is  taking  place.  The 
mathematical  statement  of  this  is  that 

V  •  n  =  0  ,  (2.3) 

where  V  •  n  is  the  normal  velocity  of  the  boundary.  However,  during 
gasification  the  temperature  at  the  surface  is  just  the  gasification 
temperature 

f ( u )  =  0  ,  (2.4) 

and  the  surface  moves  with  the  normal  velocity 

-XV  •  n  =  F  -  q—  f(u)  (2.5) 

3n 

where  X  is  the  increase  in  enthalpy  per  unit  volume  attained  by  the  solid 


unon  its  being  converted  to  gas 


There  is  no  loss  of  generality  in  choosing  u  for  the  solid  to  be  0  at 
the  gasification  temperature  0.  Then  for  the  gas  u  =  X.  We  tnav  thus  extend 
the  function  f(u),  so  far  defined  for  u  <  0,  to  the  interval  u  6  [0,X]  by 

f ( u )  =0,  0  <  u  <  X  .  (2.6) 

When  the  boundary  conditions  (2.4)  and  (2.5)  apply,  the  problem  looks  very 
much  like  a  one-phase  Stefan  problem  with  sources,  and  it  is  known  that  a  more 
concise  way  of  writing  (2.1),  (2.4),  and  (2.5)  is  in  the  form  of  a 
"conservation"  law  [1] 

ufc  =  Af(u)  +  o(x,t)  ,  (2.7) 

where  (2.7)  now  holds  over  all  space  and  f  has  been  extended  by  (2.6).  for 
our  problem,  a  has  the  form 

«*'*>  =  F53(u<0}  *  (2‘8) 

The  form  (2.7)  can  also  be  used  as  the  basis  for  an  efficient  numerical 

solution  of  N-dimensional  problems  which  avoids  the  need  for  following  the 

moving  boundary.  (The  matter  of  the  explicit  appearance  of  the  boundary  in 

the  source  term  can  be  treated  by  a  method  to  be  described  in  the  sequel.) 

There  are,  however,  some  important  limitations  to  the  applicability  of 
the  formulation  (2.7).  The  "one-phase"  Stefan  problem  is  really  a  two-Dhase 
problem  for  which  heat  conduction  takes  place  in  only  one  phase.  The  non¬ 
conducting  phase  (in  this  case,  the  "upper”  one)  acts  passively  as  long  as  the 
reoion  it  occupies  is  increasing  with  time,  but  when  the  reqion  occupied  by 
the  conducting  phase  increases  anywhere,  the  non-conducting  phase  acts  as  a 
reservoir  of  either  positive  or  negative  (in  this  case,  positive)  enerqy.  In 
the  one-phase  Stefan  problem  under  the  influence  of  externa1  f-n  free 

boundary  can  move  in  either  way.  In  the  gasification  problem,  on  the  other 
hand,  the  uae,  once  it  blows  awav,  has  no  further  influence  on  the  solid  and 
in  particular  the  free  boundary  can  move  only  one  way,  since  the  gas  cannot 


resolidify.  Thus,  although  the  formulation  (2.7),  in  which  the  equation  is 
extended  to  the  non-conducting  phase,  is  quite  natural  for  the  one-phase 
Stefan  problem,  in  our  problem  the  notion  of  a  coexisting  "gas  phase"  is  a 
fiction. 

Nevertheless,  the  fact  that  at  times  the  gasification  problem  can  be  put 
in  the  form  (2.7)  and  the  numerical  simplifications  brought  about  by  this 
possibility  serve  as  an  inducement  for  us  to  find  an  appropriate  version  of 
the  problem  in  the  spirit  of  (2.7).  We  have  observed  that  the  process  of 
gasification  is  irreversible,  so  that  V  •  n  <  0.  Thus,  we  can  summarize  the 
boundary  conditions  (2.2)-(2.5)  by  noting  that  (2.5)  holds  in  all  cases,  but 
that  on  portions  of  the  boundary,  for  certain  times,  we  have 

V  •  n  <  0,  f(u)  =  0  ,  (2.9a) 

whereas  on  the  remaining  portions  of  the  boundary  or  for  the  remaining  times, 

V  •  n  =  0,  f ( u )  <  0  .  (2.9b) 

From  V  •  n  <  0,  we  get 

q-  f ( u )  <  F  ,  (2.10) 

an 

3f 

which  looks  mathematically  something  like  an  "obstacle"  condition  on  — .  We 

on 

note  that  the  conditions  (2.9a)  are  of  Dirichlet  type  for  f  and  the 
conditions  (2.9b)  are  of  Neumann  type  for  f.  Thus,  generally  we  may  expect  a 
switching  back  and  forth  between  Neumann  and  Dirichlet  boundary  conditions, 
and  an  alternate  stopping  and  starting  of  the  boundary  motion. 

Immediately  we  discern  a  lack  of  regularity  for  the  solution  of  the 
gasification  problem  as  opposed  to  the  Stefan  solution.  por,  in  the  case  of 
the  latter,  the  temperature  f ( u)  has  been  seen  to  be  continuous  [31 .  Put  if 
we  extend  f(u)  to  the  region  where  there  is  no  solid  (and  thus  \i  2  \)  by 
the  convention  f(u)  ~  o  there,  we  see  that  in  the  case  of  the  Neumann 
conditions  (2.9b),  f(u)  is  general Iv  discontinuous.  (In  fact,  this  is  true 


no  matter  what  value  we  choose  for  f(X).)  In  addition,  because  the  solid, 
upon  gasification,  is  replaced  by  a  vacuum,  which  is  a  perfect  thermal 
insulator,  the  gasification  problem  possesses  instabilities  with  respect  to 
the  initial  and  boundary  data  which  do  not  occur  for  the  Stefan  problem.  For 

example,  by  changing  the  initial  enthalpy  to  X  on  a  set  of  arbitrarily  small 

measure,  we  can  partition  the  solid  into  a  collection  of  non-interacting 

regions,  on  the  houndary  of  each  of  which  the  boundary  conditions  (2.9b) 

obtain.  Then,  invoking  the  lack  of  continuity  of  f( u)  under  the  influence 
of  these  boundary  conditions  which  was  just  observed,  we  may  construct 
solutions  for  which  there  is  no  uniform  continuity  of  f(u)  as  the  measure  of 
the  set  on  which  u  has  been  set  to  X  goes  to  0,  and  for  which  thermal 
contact  between  regions  on  different  sides  of  the  discontinuities  of  f  has 
been  prohibited.  In  like  manner,  by  shifting  very  slightly  the  positions  of 
the  sources  of  external  radiant  energy  and  the  directions  in  which  they 
radiate,  one  may  gasify  all  the  material  along  certain  rays  and  thermally 
insulate  different  parts  of  the  solid  from  one  another,  thereby  profoundly 
affecting  the  solution.  And  by  allowing  the  radiation  to  arrive  in  sharp 
pulses  rather  than  uniformly  distributed  over  a  small  interval  of  time,  one 
may  increase  the  gasification  of  the  material  at  the  surface  while  minimizing 
the  conduction  of  heat  to  the  interior. 

In  our  picture  the  gas  has  enthalpy  u  =  X  and  is  transparent  to 
radiation,  and  hence  there  is  no  mechanism  wherebv  u  can  exceed  X.  For  the 
St°fan  problem,  one  can  give  meaning  to  a  value  of  u  in  (0,X):  it  is  that 
a  volume  fraction 


u 


+ 


\ 


v  ( ') )  =  -  max  (  u  ,  h  ) 
A 


(2.11) 


at  a  point  has  enthalpy  X  ami  is  gas,  while  the  remainder  is  solid  with 
enthalpy  0.  Since  values  of  u  in  ( 0 , X )  arise  quite  naturally  when  one 
solves  the  Stefan  problem  with  sources,  (2.7),  we  may  expect  a  similar 
situation  to  develoo  when  we  try  to  cast  the  gasification  problem  in  like 
form.  However,  the  values  permitted  the  enthalpy  at  the  solid  boundary  are 
now  not  restricted  to  the  set  {X, 0},  but  may  belong  to  the  set  {X,u  <  0}. 
Accordingly,  if  we  denote  the  volume  fraction  at  each  point  that  corresponds 
to  enthalpy  A  by  x»  so  that  ( 1  -  x)  denotes  the  volume  fraction  of  solid 
with  enthalpy  u  <  0,  we  get  for  the  total  enthalpy  u  of  the  combination, 

u=u(1-x>+xX»  (2.12a) 

or 


u  -  u 

x  = - 1 

X  -  u 


>  X(u) 


(2.12b) 


Here  we  have  proceeded  as  if  the  solid  material  had  a  common  enthalpy  u  when 

X  <  1.  We  could  envisage  situations  in  which,  in  the  region  where 

0  <  x  <  we  bad  volume  fractions  on  of  solid  material  with  enthalpies 

u.  <0,  each  set  thermally  insulated  from  the  others.  Then  we  would  have 

l  a.  =  1  -  x,  l  a.  u  =  (1  -  x)u  .  (2.13) 

i  i 


However,  if  we  consider  the  assignment  of  initial  and  boundary  conditions  and 
the  gasification  process  itself  as  phenomena  which  are  essentially  subject  to 
mul t i- d imens iona 1  stochastic  spatial  fluctuations,  albeit  minute  ones,  it 
annears  that  the  comnartmentalization  of  the  solid  mater  ill  on  a  microscopic 
level  into  such  components  thermally  isolated  from  each  other  would  occur  with 
probability  0,  and  we  shall  consider  the  solid  material  to  have  a  sinale- 
••alued  enthalpy  u  and  si nole-va lued  temperature  flu). 
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From  the  foregoing  considerations,  we  may  derive  equations  for  the  time 


evolution  of  x*  If  gasification  is  taking  place  we  have  u  =  0  and 


X  =  X<u> '  \  *  < < x<u> )t>+ 


(2.14a) 


But  if  gasification  is  not  takinq  place,  we  have  Neumann  conditions  at  the 
solid  boundary,  u  <  0,  and 

X  >  X<u>  <•  Xt  =  0  •  (2’  14b) 

If  x(x,0)  =  y(u(x,0))  Vx,  we  can  use  (2.14a,b)  to  find  y(x,t)  immediately 
in  terms  of  u: 

y(x,t)  =  sup  •  (2.15) 

0  <t  *  <t 

As  regards  the  absorption  of  energy  at  the  solid  surface,  as  given  by 
(2.9),  we  treat  the  gas  as  if  it  were  completely  transparent  to  the 
radiation.  The  absorption  of  radiation  in  a  region  where  gas  and  solid  are 
interspersed  on  a  microscopic  level  should  then  be  proportional,  at  each 
point,  to  1  -  x*  Accordingly,  if  0  denotes  the  location  of  the  source  of 
radiation,  and  we  assume  the  source  to  be  isotropic,  a  relevant  quantity  at 
any  point  P  is  the  following  integral: 

0 

d(P,0>  =  /  (1  -  x(Q))ds  ,  (2.16) 

P 

where  the  integral  is  taken  over  the  straight  line  connecting  P  to  0.  We 
consider  the  soLid  to  be  completely  opaque  to  the  ra  nation,  and  leads  to 

the  requirement  that  all  the  radiant  enerov  be  deposited  at  points  p  where 
d(P,f>)  is  as  small  as  possible,  subject  only  to  the  constraint  that 

u(x,t)  <  1  Vx , t  .  (2.17) 
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We  have  presented  the  physical  considerations  wh ;  hiv<>  m.  •••  «  :  ir 

treatment  of  the  gasification  problem.  In  the  next  section  we  «  i  i  i  nreserr  an 
algorithm  to  solve  a  time-discretized  version  of  the  oroblem.  ;s  al  mu;1  s 
will  be  of  a  sort  that  has  been  introduced  to  solve  problems  of  the  tyne  (2.7/ 
[1].  We  should  note,  however,  that  although  we  may  resort  at  times  in  the 
sequel  to  pseudo-physical  language  to  interpret  the  algorithm  which  we 
develop,  the  algorithm  to  be  presented  will  have  been  crafted  to  solve 
precisely  the  problem  whose  physical  and  mathematical  countenances  have  been 
so  far  explained.  The  algorithm  should  not  necessarily  be  construed  as 
describing  also  any  physically  realistic  problems  of  a  more  general  nature. 
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3.  Time-Discretized  Formulation 


If  the  Lipschitz  coefficient  of  f(u),  as  given  in  (2.1),  is  1/f?,  a 
suitable  numerical  scheme  to  solve  (2.7)  with  0=0  is,  with  t  a  time  step 
and  u(nr)  approximated  by  un  [1]  , 

un+1  =  un  -  6f(un)  +  S(|)(6f(un))  ,  (3.1) 

D 

where 

S(h)  H  ehA  (3.2) 

is  simply  the  linear  semi-group  generated  by  the  operator  A.  In  the  problem 
we  treat  here,  we  mav  regard  S(h)  as  Green's  function  for  the  heat  equation 
in  RN.  It  will  be  convenient  below  to  write  out  S(h)  in  this  fashion 
explicitly : 

( S(h)v){x)  =  f  S(h;x,x* )v(x' )dx’  .  (3.3) 

N 

R 

A  pseudo-physical  interpretation  of  (3.1)  is  that  the  equation 

uv  =  ~  A(ff?)  (3.4) 

p 

is  an  equation  giving  the  evolution  of  u  in  terms  of  a  diffusion  of  (3f(u), 
and  that  to  find  un+ 1  we  take  un,  subtract  off  f?f(un),  which  is  to  be 
diffused,  and  then  add  back  S(-j)  (  fjf(un)  )  ,  the  result  of  the  diffusion. 

If  we  are  dealing  with  a  problem  for  which  (2.T4a)  applies  Vx,t  such 
that  0  <  x  <  1,  then  in  fact  (3.1),  supplemented  by  a  term  to  represent  the 
effect  of  the  external  heat  sources,  will  satisfactorily  describe  the  time 
evolution.  However,  in  the  more  general  case  in  which  we  may  also  have 
X  >  <(u<  and  Neumann  conditions  at  the  boundary,  this  algorithm  clearly  will 
not  suffice. 

In  the  first  place,  according  to  the  pseudo-pbysic3l  type  of  reasoning 
introduced  above,  the  true  measure  of  the  "thermal  enerqy"  of  the  solid,  per 
unit  volume  of  solid,  should  be  ,3f(u)  and  not  Bf(u),  and  the  total  amount 


r 


of  such  thermal  enercry,  oer  unit  volume,  is  Sf(u)(1  -  y).  Here,  in  accord 
with  (2.1 2a) , 


u 


u  2  XX 
1  -  X 


(3.5) 


when  y  <  1.  Accordingly,  we  would  subtract  out  0f(u)(1  -  x) »  then  add  back 

S(-j)  [  (3f  (u)  ( 1  -  v)].  However,  in  diffusing  the  thermal  energy  (3f(u)(1  -  y)  < 

P 

we  have  transmitted  energy  to  regions  where  a  fraction  y  of  the  volume  acts 

as  a  perfect  insulator.  Hence  this  portion  of  the  diffused  thermal  energy 

should  be  subtracted  out  and  returned  to  the  point  from  whence  it  came.  This 

means  subtracting  out  yS(-~)  [  0f  (  u)  ( 1  -  y)]  and  then  adding 

f  S(-~; x' ,x) y(x' ) Bf (u( x) ) ( 1  -  y(x))dx*.  This  final  subtraction  of  a  portion 
R  B 

of  the  diffused  thermal  energy  and  return  to  its  source  is  the  part  of  the 
algorithm  designed  to  conform  to  the  insulating  properties  of  the  vacuum  and 
to  yield  the  solution  of  a  Neumann  problem  at  the  appropriate  places  on  the 
boundary.  It  is  a  direct  generalization  of  the  method  of  "images"  used  to 
solve  Neumann  problems  at  planar  boundaries  for  the  classical  equations  of 
mathematical  physics.  Putting  these  arguments  together,  for  the  case  when 
Green's  function  is  symmetric,  as  it  is  here, 

S ( h ; x , x ’ )  =  S ( h ; x ' ,x)  ,  ( 3.6  ) 

we  are  led  to  construct  the  function 


un  -  Bf(un)(1-xn)  +  3(-|)f,Bf(un)(1-xn)] 


(3.7; 


-n 


o  ^ 

If  there  were  no  external  energy  sources,  we  would  set  u  equal  to 

11°.  In  that  case,  on  npnroximatinq  S(h)  by  the  first  two  toms  of  its 
fornul  Taylor  expansion. 
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S(h)  -*■  1  +  hA 


(3.8) 


we  would  get  from  (3.7),  as  t  -*■  0, 

ut  =  (1-x)A((1-x)f(u))  -  d-x)f(u)  A  ( 1  -  x )  •  (3.9  ) 

In  fact,  (3.9)  may  be  given  a  pseudo-physical  interpretation  by  noting  that  it 
can  be  derived,  in  the  limit  as  the  "mean  free  path"  p  *  0,  from  the 
Boltzmann- like  equation 

Up  =  ~  { /  g(~-X---)f  (u(  x'  ) )  ( 1-x(x) )  ( 1-x<x'  ))dx' 

rN  p 

(3.10) 

-  f(u(x))(1-x(x))  /  g^-7-'-1  )(1-x<x'))dx'  }, 

VT  P 


o 

where  /  g(  r)  dr  =  1 ,  g  >  0,  and  Ap  =0(1)  as  p  ■+•  0.  Formally,  (3.9)  may 
0 

also  be  written 

ut  =  7  •  ((1  -  x>2Vf(S>>  .  (3.11  ) 

However,  when  f(u)  and  x  are  discontinuous,  as  may  be  the  case  at  a  sharp 
solid-gas  interface  in  the  circumstance  of  Neumann  conditions  at  the  boundary, 
the  meaning  of  (3.11)  in  a  distributional  sense  is  not  clear. 

In  the  general  case  we  have  to  add  the  effect  of  the  external  energy 
sources  to  3°  in  order  to  get  un+^.  Let  us  suppose  that  the  incident 
energy  F  in  (2.8)  has  the  form 


F(  x ,  t ) 


F(t)  |n  «  (x(Q,>  -  x)  I 

Vi  | x(0  )  -  x|N 


(3. 12) 


where  x(0)  is  the  location  of  the  source  and  =  |3B^|,  the  unit 

hall  in  N  dimensions. 
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The  energy  deposited  between  times  nt  and  (n  +  1)t  may  reasonably  be 


approximated  by  T0(x,nr)  when  0  is  continuous  with  respect  to  t.  When 
0  has  the  form  (2.8),  this  quantity  has  a  singularity  at  the  solid 
boundary.  However,  during  the  time  t  that  energy  is  deposited,  it  is  also 
being  conducted  into  the  interior  in  the  manner  described  by  equation  (2.1). 
Accordingly,  if  fi  is  a  local  value  of  f'(u)  in  the  solid  and  fi  >  0, 


1>V2 


the  energy  will  be  spread  out  over  a  distance  0( ( i£n  — )  2  )  after  the  time 
interval.  It  will  be  important  for  us  to  take  this  fact  into  account  when 
depositing  the  incident  radiant  energy.  To  spread  the  energy  out  over  a 

1  Vo 

distance  larger  than  0( ( tin  — )  2  )  would  be  needlessly  to  increase  the  width 
of  the  free  boundary  (that  is,  the  distance  over  which  x  changes  from  "near 
0"  to  "near  1")  generated  by  the  algorithm.  Experience  with  the  Stefan 
problem  leads  us  to  expect  an  approximate  free  boundary  of  width 


n+1 


0( ( x£n  — )  z  ) ,  anyway  [2].  Thus,  we  calculate  u  from  the  equation 


n+1 


(  x)  =  0°  (  x)  + 


/t  F(  m)  ( 1-y  (x)  ) 

|x(°)-x|N  1 


(3.  13a) 


0  <  d  ( x)  <  /t 
dn(x)  >  / t 


where,  from  (2.16), 


0 

dn(  x)  =  /  (1  -  xn(Q>><ls  • 

X 


(3. 13b) 


r’nee  nn+1  ^as  been  found , 

semi-discrete  version  of  (2.15): 

-n+1 

:< 

I f  we  refer  to  (2.12a)  and  (3.6), 


X  is  determined  in  accordance  with  a 

max ( x” » u°+1 ) )  .  (3.14) 

we  find  from  (3.7)  that 
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and  thus 


rfn  - 


<  0 


1  - 


xifi  <  Xn  • 


(3.15) 


Accordingly,  if  there  are  no  external  energy  sources  and  un+1  =  Cin,  we  get 
from  (3.14)  that  yn+ 1  =  yn.  This  simply  says  that  no  material  is  gasified  in 
the  absence  of  external  sources. 

We  have  only  to  check  that  the  algorithm  (3.7),  (3.13),  (3.14)  guarantees 

that,  if  0  <  x°  *  1 /  then  the  following  two  conditions  are  satisfied: 
n+ 1 

(i)  0  <  y  <1  and  (ii)  the  time-discretized  version  of  (2.17)  holds,  or 

n  _i_  4 

u  <  X.  These  two  conditions  are  actually  equivalent,  by  (2.11),  (3.14), 
and  the  assumption  0  <  y11  *  By  direct  calculation  from  (2.12a),  (3.7), 
and  (3.13),  we  get 


u 


n+1 


<  xnx  + 


( 1  -  xn)  /t  F(  nr) 

Vl|x(0)  -  x|N-1 


(3. 16) 


(3.16)  will  imply  that  un+1  <  X  if  we  impose  the  following  stability 
condition  on  the  size  of  the  time  step  t: 

F(nx)A  <  Xo^_1  (dist(0 ,  {x|u(  x,0 )  <  X})]N  1  .  (3.17) 

The  sufficiency  of  (3.17)  follows  from  the  observation  that,  if  y11  (x)  <  1, 
then  y^(x)  <  1  and  |x(0)  -  x|  >  dist(0,  {x|u(  x,0 )  <  X}).  In  the 
calculations  to  be  described  in  the  next  section,  we  impose  the  stability 
condition  (3.17)  on  t. 
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4.  Numerical  Quadrature  and  Computational  Results 


The  computations  we  have  actually  performed  have  been  one-dimensional 
calculations  for  the  case 

X  =  1,  f(u)  =  min(u,0)  .  (4.1) 
We  see  from  (4.1)  that  we  can  choose  g  =  1.  The  algorithm  described  in  the 
last  section  can  now  be  written 


un+  =  un+(1-xn/[(S(T)-1)(1-xn)f(un)+f(un)(S(i)-1)xn] 

[  /t  (1  -  xn>  o  <  dn  <  A  (4.2a) 

+ 

0  dn  >  /t 


-n+1  ,-n  ,  n+1  ,  .  ... 

X  =  max( x  #X<U  > )  •  (4.2b) 

Since  the  expression  (3.5)  for  u  becomes  indeterminate  when  x  *  we  have 

used  the  following  determinaton  of  (1  -  that  is  well-defined  for  all 

values  of  x  e  [0,1]: 

(1  -  x>f(u>  =  min(u  -  x»°>  •  (4.3) 

We  have  broken  the  spatial  region  into  a  finite  number  of  cells  1^  of 
width  Ax^,  1  <  i  <  I.  We  have  imposed  the  conditions  x=0,  f(u)=-1  at 
the  left  and  X  =  ^  at  the  right,  and  we  have  considered  the  source  to  be  at 
the  right  of  the  mesh.  We  denote  the  characteristic  function  of  1^  by  x-  : 

(1  x  e  I. 

1  (4.4) 

0  xei.  , 

we  aoproximate  the  functions  un(x)  and  x^*)  by  the  piecewise-constant 


functions 


(4.5a) 


n  r  n 

U  ~  l  Vi  ' 

1=1 


-n  r  -n 

X  ~  L  X^X, 
i=  1 


(4.5b) 


u"  =  (Ax.  )uj 


and  we  denote  the  total  energy  in  by  u9: 

(4.6) 

For  the  operator  S(x),  we  have  used  an  explicit  finite-difference 
scheme.  Then  the  operator  can  be  represented  by  its  effect  on  a 
characteristic  function: 


Ax . 

S(T)X.  =  T—  C 


.  0  +_^i_  + 

i  Ax  “i*i-1  Cixi  Ax  Cixi+1 
i-1  1+1 


(4.7) 


Suitable  coefficients  cT,  c9,  ct  can  easily  be  calculated,  and  we  note  only 
the  result  of  such  a  calculation  [6,  p.  29] : 


-1 


c.  =  3x[Ax._  1/  <&*i_1+&Xi+Axi+1)1  • 


(4.8a) 


c.  =  3x[  Ax  i.  (Ax.  +  Ax,  +  Ax .  )  ] 

i  1+  V2  i-1  i  1+1 


-1 


(4.8b) 


Here  we  have  set 


and 


°  ,  -  + 

c.  =  1  -  c.  -  c.  . 

i  i  i 


ci+  V,  =  2  Uxi  +  &Xi+1 


(4.8c) 


(4.9a) 


Ax^  =  Ax  1 ,  Ax  =  Ax  .  ( 4 .  Qh) 

We  recall  that  the  operator  S( x)  is  symmetric,  as  stated  in  (3.6). 
However,  the  finite-difference  analoq  of  S ( x)  given  above  is  not  necessarily 
symmetric  if  the  mesh  widths  Ax^  vary  with  i.  This  lack  of  symmetry  could 
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lead  to  a  departure  from  strict  energy  conservation  if  one  were  to  use  the 
algorithm  (4.2)  unthinkingly.  Such  a  departure  can  be  avoided  if  we  return  to 
the  argument  given  in  the  paragraph  immediately  preceding  (3.6).  What  we 
should  do  is  use  the  finite-difference  representation  of  S( r)  given  above 
when  S(t)  operates  on  (1  -  xn)f(un),  and  use  the  transpose  of  this 
representation  when  S(t)  ooerates  on  yn.  The  transpose  of  the  finite- 
difference  representation  is  given  by 

n 

(4.10) 

The  explicit  nature  of  our  representation  of  S(t)  forces  on  us  an 
additional  stability  condition  for  t.  A  sufficient  condition  is  that  t 
satisfy 

1  2 

T  <  —  mm(  (  Ax.  )  )  .  (4.11) 

2  .  l 

l 


sT<T)Xi  =  ViXi-!  +  +  Wi+i 


The  finite-difference  treatment  of  the  energy  deposition  from  outside 
apportions  to  the  cell  1^  the  amount  of  energy  V™,  given  by 


V  =  min 


-,(  1-xt)Axt/t 


■)  . 


(4.12a) 


..n  v  ,.n  ..  ~n- ■  r  nnti'n 

vi  ~  lminl  2 —  T~  /•  Vi'(1"Xi)&JVT — - — )) 

j=i+1  J  + 


1  ,,  -n.  ...  fl  ft  nx) 


1  <  i  <  I  -  1 


(4.12b) 


Denoting 

y”  H  min(u^  -  x",0)  (4.13) 

and  putting  the  expressions  above  together,  we  find,  for  2  4  i  <  T  -  1 , 
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uin+1  -  ui +  +  Yi-1ci-1Axi-13 


^[(l-x^lcj  +  (1"?-1)CI1  +  Vi  * 


(4.14) 


i 


Special  expressions  are  needed  for  i  =  1  and  i  =  I.  Since  we  are 
imposing  the  conditions  x  =  ^  and  f(u)  =  -1  at  the  left  and  X  ~  ^  at  '-'le 
right,  we  set 

—  r%  n  rs 

(4.15) 


\  '  °'  '«?♦> '  ’•  • 0 


Thus  we  get  for  i  =  1 

uf1  =  u"  +  ( 1-x”)  (Y2C2Ax2-C0  ^1 )  “  Y1  ^1  t(1-X2)Cl+C13  +  V1  '  (4.16a) 


where 


c„  = 


3t 


0  Ax1(2Axl  +  Ax2)  ' 


(4. 16b) 


and  for  i  *  I 
n+1 


-n  n  + 


and 


°i  -  ui  *  "-".'V.'i-ih-,  -  ’W’-W'i  * v! 

The  computational  loop  is  completed  by  setting 

n+1  n+1  , , 

Ui  "  Ui  /Axi 


-n+1  ,  n+1  .  -n 

x.  =  max(  ( u,  )+,x^>  • 


(4.16c) 


(4. 17a) 


(4. 17b) 


"1  i  + 

In  the  more  general  one-dimensional  case,  we  get  a  finite-difference 

relation  of  the  sort 
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un+1 

1 


=  Un  -  R(1-Y?)f(u?)Ax. 
1  A1  1  1 


-n 


+  {^xi+ici+1^(ui+1 ) ( 1-Xi+i ^ 


+  Ax,c°f(un)  ( 1~xn)  } 
111  A1 

+  f(un)(1-vn)Ax.(c+x"  , 
l  Hi  iH+1 


0-n,  — n  v  ,  ,,n 
+  c.  x.  +c.  x.  .  )  + 
i  H  i  H- 1  l 


(4.18> 


where  and  cT  are  as  given  in  (4.8a,b),  and 

ci  =  6  ~  ci  '  ci  *  (4.19) 

A  sufficient  condition  for  stability  of  the  finite-difference  scheme  is 

t  <  r  6  min{ Ax. )2  .  (4.20) 

i  1 

is  given  in  terms  of  u?  and  ^  by  means  of  (3.5).  For  computational 
reasons,  for  a  sufficiently  small  e,  one  may  wish  to  calculate  u™  from 


-n 
u . 
1 


n 

u. 

i 


x*i 


-n 

X, 


1 


e 


(4.21) 


0  1  -  x"  <  e 

We  performed  a  calculation  with  I  =  120,  Ax^  =  .01  Vi,  x  =  .00005, 
u?  =  -1,  1  <  i  <  100,  =  1,  101  <  i  <  120,  and  x^  =  chose 

F(t)  =6  independent  of  t. 

Letting  the  intersection  of  *100  and  be  located  at  x  =  0  and 

denoting  by  the  "exact  solution"  the  solution  of  the  gasification  oroblem 
nosed  on  R  with  the  solid  region  initially  occupying  R_  and 
u(x,0)  =  -1  V*  fi  R_,  we  find  that  initially  the  exact  solution  satisfies 
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Neumann  conditions  at  x  =  0,  and  that  u(0,t)  is  given  by  -1  +  6/t/ir, 
until  the  time  t^  =  n/36,  when  gasification  commences  and  the  front  begins 
to  move.  Further  analysis  of  the  exact  solution  indicates  that  for  t  >  tQ 

but  t  -  tQ  small,  the  position  of  the  front  is  at 

54  2 

Ht)  =  -  —  (t  -  tQ)  . 

IT 

If  the  solution  of  the  gasification  problem  has  a  self-similar  behavior  as 
t  >  »,  it  will  be  of  the  form 

|  x+  J  Ct-t  ) 

u  +  -1  +  e  for  x  <  -  j  (t  -  t,)  . 

The  problem  was  run  to  time  t  =  .3.  The  following  table  gives  the  computed 
position  £(t)  of  the  free  boundary,  defined  to  be  the  point  where  y  =  .5 
and  determined  by  linearly  interpolating  values  of  x  between  adjacent  cells 
such  that  x  “  *5  changes  sign  from  one  to  the  other.  In  all  cases  the  front 
was  never  spread  out  over  more  than  one  or  two  computational  cells. 


TABLE  4. 1 


Computed  Position  E(t)  of  Free  Boundary 


t 


lit) 


t 


Z(t) 


09 

0 

10 

-.005 

11 

-.010 

12 

-.019 

13 

-.028 

14 

-.037 

15 

-.046 

16 

-.055 

17 

-.064 

18 

-.076 

19 

-.087 

.20 

.21 

.22 

.23 

.24 

.25 

.26 

.27 

.28 

.29 

.30 


-.098 
-.  109 
-.  1  19 
-.  130 
-.  142 
-.  154 
-.  166 
-.  177 
-.  189 
-.200 
-.212 
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In  order  to  test  that  the  algorithm  described  here  does  indeed  predict 

the  cessation  of  gasification,  cooling  off  of  the  solid  surface,  and 

resumption  of  Neumann  conditions  at  the  boundary  under  the  appropriate 

circumstances,  for  m  >  .3  we  set  u?,  =  0.  As  expected,  the  motion  of  the 

6 1 

free  boundary  soon  ceased,  and  the  steady  state  value  of  the  position  was 
computed  as  E  =  -.216.  The  solution  between  i  =  61  and  the  free  boundary 
reached  a  steady  state  u  with  the  slope  ux  =  3,  as  it  should,  from  (2.2) 
and  (3.12). 
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5.  Renarks  on  Monotonicity 

We  would  like  to  conclude  this  paper  with  some  suggestions  regarding  the 
sorts  of  monotonicity  properties  we  might  expect  for  the  solution.  It  is  well 
known,  of  course,  that  the  solution  of  the  Stefan  problem  depends 
monotonically  on  the  initial  and  boundary  data,  and  also  on  the  external 
sources.  In  the  case  when  y  =  y(u)  Vx,t,  then  we  also  have  such  monotone 
dependence.  For,  we  are  then  dealing  with  the  situation  in  which  the 
gasification  front  is  moving  and  f(u)  =  0  at  the  free  boundary.  In  this 
event,  (3.11)  is  meaningful  in  a  distributional  sense  and  may  be  rewritten  as 

ufc  =  A0( u)  (5.1a) 

for  an  appropriate  function  9: 

9'  ( u)  =  (1  -  X(u))2f'(u)  .  (5.1b) 

But  (5.1a)  is  in  the  same  form  as  the  Stefan  problem,  and  the  monotonicity 
follows  in  the  usual  way. 

In  general,  we  can  show  that  this  type  of  monotonicity  does  not  hold  for 
the  gasification  problem.  Consider  a  problem  with  f(u)  =  min(u,0)  and 
\  =  1,  and  with  the  boundary  condition  (f(u))x  -  0  imposed  at  x  =  0. 
Suppose  the  initial  data  are  u(x,0)  =  -1  for  0  <  x  <  1,  u(x,O)=0  for 
1  <  x  <  2,  and  u(x,0)  =  1  for  x  >  2.  If  we  do  not  irradiate  the  boundary 
of  the  solid,  at  x  =  2,  with  energy,  the  problem  will  evolve  as  a  diffusion 
process  in  0  <  x  <  2  with  homogeneous  Neumann  data  at  x  =  2.  Eventually 
the  steady  state  u(x)  ~  for  0  <  x  <  2,  u( x)  =  1  for  x  >  2  will  be 

achieved.  On  the  other  hand,  suppose  we  irradiate  the  solid  with  =  — 

2  o 

for  0  <  t  <  6  and  F(t)  =0  for  t  >  6.  As  we  let  6+0,  we  find  that 
any  diffusion  of  thermal  energy  in  the  region  occupied  by  the  solid  in  the 
time  interval  0  <  t  <  6  becomes  negligible.  Then,  in  the  limit  as  6+0, 
we  have  for  u(x,6)  the  following  values:  u(x,6)  =  -1  for  0  <  x  <  1, 
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u(x,6)  =  1  for  x  >  1.  This  is  also  the  steady  state  achieved  for  the 
problem.  Thus,  by  irradiating  the  solid  surface  with  energy,  we  have  actually 
reduced  the  steady  state  values  of  u  in  portions  of  the  interior. 

A  tempting  conjecture  to  make  is  that,  if  we  have  two  solutions 
uf(x),  x>)  and  UjCx),  ,£(x)  of  the  gasification  algorithm,  subjected  to 
the  same  boundary  conditions  and  external  energy  source  located  at  a  point 
with  x-coordinate  equal  to  X,  then  if  u^{x,y)  >  u^XjyJVXjy  and 
X"(x,y)  =  x;<*,y>Vx.y,  it  will  follow,  for  any  p  >  0,  that 

X+J,  X+f 

/  /  u"+P(x,y)dxdy  >  f  f  u"+P(x,y)dxdy 

rN-1  X-C  rN-1  x-c 

and 

X+£  X+£ 

/  /  X?  P(x,y) dxdy  >  f  f  Xo  P(x,y)dxdy 

rn-i  X-C  rn-i  x-5 

for  all  £  e  R+.  Here  y  merely  labels  the  N  -  1  coordinates  of  a  point  in 
a  hyperplane  orthogonal  to  the  x-axis.  If  the  conjecture  is  true,  we  will  get 
monotonicity  properties  for  Radon  transforms  of  the  energy  distributions. 
However,  at  this  point  we  do  not  have  any  solid  indications  regarding  the 
accuracy  of  the  conjecture. 
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